Where post-Newtonian and numerical-relativity waveforms meet 
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We analyze numerical-relativity (NR) waveforms that cover nine orbits (18 gravitational-wave 
cycles) before merger of an equal-mass system with low eccentricity, with numerical uncertainties 
of 0.25 radians in the phase and less than 2% in the amplitude; such accuracy allows a direct 
comparison with post-Newtonian (PN) waveforms. We focus on one of the PN approximants that 
has been proposed for use in gravitational-wave data analysis, the restricted 3.5PN "TaylorTl" 
waveforms, and compare these with a section of the numerical waveform from the second to the 
eighth orbit, which is about one and a half orbits before merger. This corresponds to a gravitational- 
wave frequency range of Mlo = 0.0455 to 0.1. Depending on the method of matching PN and NR 
waveforms, the accumulated phase disagreement over this frequency range can be within numerical 
uncertainty. Similar results are found in comparisons with an alternative PN approximant, 3PN 
"TaylorT3". The amplitude disagreement, on the other hand, is around 6%, but roughly constant 
for all 13 cycles that are compared, suggesting that only 4.5 orbits need be simulated to match PN 
and NR waves with the same accuracy as is possible with nine orbits. If, however, we model the 
amplitude up to 2.5PN order, the amplitude disagreement is roughly within numerical uncertainty 
up to about 11 cycles before merger. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 98.80.Jk 



o 
cr. 



> 

in 
o 

m 

o 
i> 
o 



X 



I. INTRODUCTION 

The current generation of interferometric 
gravitational-wave detectors [H, 0, H[ have reached 
design sensitivity, and have recently completed taking 
data in the S5 science run. Signals from coalescing 
black-hole binaries will be among the strongest that one 
hopes to find in the detector data, and data analysts 
are searching for them by performing matched filtering 
against template banks of theoretical waveforms. At 
present data analysts use, or are preparing to use, 
waveforms calculated by post-Newtonian (PN) meth- 
ods, and in particular the standard Taylor-expanded, 
effective-one-body (EOB) and BCV [4j] waveforms 
implemented in the LSC Algorithms Library (LAL) [j|, 
although current GW searches do not go beyond second 
PN order (2PN). The PN waveforms are expected to 
be reasonably accurate during the slow inspiral of the 
binaries, but it is not clear how well they can model 
the merger phase. Ultimately the PN waveforms will be 
connected to waveforms from fully general relativistic 
numerical simulations, which will also model the last 
orbits, merger and ringdown. 

In the last two years breakthroughs in numerical rel- 
ativity [H, 0, IH have completed the work of providing 
the techniques to generate the necessary numerical (NR) 
waveforms. The nonspinning equal-mass case in partic- 
ular has been studied in great detail [1, @, H, H, OH [HI, 
[H, [H, [3, [ID, [H, [TtJ and extremely accurate waveforms 
over many (> 15) cycles are now available [HI]. A first 
comparison of PN and NR equal-mass waveforms was 
made i n flBL [l9| . unequal-mass waveforms were studied 
in [2(J, l2ll. |22|. and spinning binaries in (23| , and the 
work of producing hybrid NR-PN waveforms has begun 



|21l |22j, |24|, |25j . Good agreement has been observed be- 
tween NR and PN waveforms [HI, [n| Hi]], and in par- 
ticular phase disagreements of less than one radian up 
to ~ 1.5 orbits before merger were seen in [ll|. How- 
ever, until now NR waveforms have not been accurate 
enough to allow a conclusive comparison with the PN 
wave amplitude; for example, it was pointed out in (2~l| 
that although the disagreement in the amplitude of NR 
and PN waves was about 10%, this was also the size of 
the uncertainty in the NR wave amplitude, and it was 
not possible to conclude what order of PN treatment of 
the wave amplitude gives the best agreement with fully 
general-relativistic results. 

In this work we systematically compare numerical 
equal-mass waveforms that include up to 18 cycles before 
merger with the 3.5PN "TaylorTl" and 3PN "TaylorT3" 
waveforms implemented in LAL. One could compare with 
many different varieties of PN waveform, but the Tl and 
T3 approximants are common choices that are among 
those proposed for gravitational-wave searches in detec- 
tor data, and restricting ourselves to only two approxi- 
mants keeps our analysis and the presentation of our re- 
sults relatively simple. The region of comparison includes 
13 cycles. Considering the amplitude A(t) and phase <j){t) 
of our numerical waveforms separately, we find that the 
accumulated error in <f>(t) is at most 0.25 radians over 
the frequency range of comparison. These uncertainties 
are dominated by the finite extraction radii of our wave- 
forms, not finite-difference errors. The error in the am- 
plitude A{t) is less than 2% for most of the simulation. 
We estimate the eccentricity as e < 0.0016. We therefore 
consider these waveforms to be adequately accurate for 
a detailed comparison with PN results, in particular to 
determine the disagreement between NR and PN wave 
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amplitudes. 

Numerical simulations are computationally expensive, 
and mapping the parameter space of binary mergers (in- 
cluding black holes of varying mass ratio and spins) will 
require huge computer resources. As such, we would 
prefer to simulate only a small number of orbits before 
matching with PN results. We find that a simulation of 
only 4.5 orbits has the same amplitude agreement with 
the last four cycles of the restricted 3.5PN waveform 
as the long 18-cycle simulation, and therefore suggest 
that relatively short numerical simulations are feasible 
for matching to PN inspiral waves. For an even greater 
amplitude agreement with PN theory, our results suggest 
that, using 2.5PN amplitude corrections, at least 5.5 or- 
bits (11 cycles) before merger are necessary. 

We also compare the black holes' motion with that 
calculated by integrating the PN equations of motion 
[H, H3 and find that the PN and NR orbital tracks and 
frequencies are in excellent agreement until the last three 
orbits of the binary. 

Before describing our analysis in detail, we give a 
brief summary of our numerical methods in Section |TT] 
and the procedure for generating PN waveforms in Sec- 
tion IIIII In Section IIVI we discuss the simulations we 
performed, and in particular establish the sixth-order 
convergence of our results, construct Richardson extrap- 
olated waveforms with error estimates, and extrapolate 
the finite-extraction-radius waveforms to those measured 
as R ex — * oo. We also discuss the phase errors in our 
waveforms and give a consistency check between waves 
from simulations starting at different initial separations. 
In Section [V] we directly compare the PN and numerical 
waveforms. 



TABLE I: Summary of grid setup for numerical simulations. 



Grid 


hmin 


hmax 


'max 


D8 simulation 


X„=2[5 x 56 : 5 x 112 : 6] 


M/37.3 


m/7M 


775M 


D9 simulation 


X„=2[5 x 56 : 5 x 112 : 6] 


M/37.3 


m/7M 


775M 


D10 and Dll simulations 


Xt)=2 [5 x 48 : 5 x 96 : 6] 
X„=2[5 x 56 : 5 x 112 : 6] 
X„=2[5 x 64 : 5 x 128 : 6] 


M/32.0 
M/37.3 
M/42.7 


16M 
96/7M 
12M 


776M 
775M 
774M 


D12 simulations 


X„=2[5 x 64 : 5 x 128 : 6] 
X„=2[5 x 72 : 5 x 144 : 6] 
X„=2[5 x 80 : 5 x 160 : 6] 


M/42.7 
M/48.0 
M/53.3 


12M 
32/3M 
48/5M 


774M 
773M 
773M 



TABLE II: Physical parameters for the moving-puncture sim- 
ulations: the coordinate separation, D/M, the mass param- 
eters in the puncture data construction, rrii/M, and the mo- 
menta p x /M and p y /M. The momenta are based on those de- 
scribed in [53], and produce quasi-circular inspiral with min- 
imal eccentricity. The estimated eccentricity e is also given, 
as described in the text. The punctures are placed on the y- 
axis, and for all simulations the total initial black-hole mass 
is M = 1. 



Simulation 


D/M 


rm/M 


Px/M 


Py/M (X10- 3 ) 


e 


D8 


8.0 


0.48240 


T-0.11235 


T2.0883 


0.0025 


D9 


9.0 


0.48436 


TO. 10337 


+1.4019 


0.0022 


D10 


10.0 


0.48593 


T0.096107 


T0.980376 


0.0022 


Dll 


11.0 


0.48721 


T0.090099 


T0.709412 


0.0020 


D12 


12.0 


0.48828 


T0.085035 


T0.537285 


0.0016 



II. NUMERICAL METHODS AND 
WAVEFORMS 

We performed numerical simulations with the BAM 
code [ill, , replacing fourth-order accurate derivative 
operators by sixth-order accurate spatial derivative op- 
erators in the bulk as described in [lj|. The code starts 
with black-hole binary puncture initial data [2!| [3(| gen- 
erated using a pseudo-spectral code [3l|, and evolves 
them with the \~ variant of the moving-puncture [32l |33| 
version of the BSSN 0, ilformulation of the 3+1 Ein- 
stein evolution equations [361 ] - The gravitational waves 
emitted by the binary are calculated from the Newman- 
Penrose scalar ^4, and the details of our implementation 
of this procedure are given in [l3| . 

The simulations we performed for this analysis are 
summarized in Tables U and [TTJ For the configurations 
with initial separations D = 1QM, 11M, 12M (denoted 
by "D10", "Dll" and "D12" throughout the paper), sim- 
ulations were performed at three resolutions, and final re- 
sults obtained by Richardson extrapolation, as described 
in Section[lV] For the D = 8M, 9M ("D8" , "D9") config- 
urations, which are used only for comparison at the end 



of Section [Vl only one simulation at medium resolution 
was performed. 

The physical parameters are given in Table [TTJ The 
initial momenta for low-eccentricity quasi-circular inspi- 
ral are estimated by the PN method described in [271 ]. 
We estimated the eccentricity from the frequency io v of 
the puncture motion, as we did previously for D = 11M 
simulations in [13] , and as also used in [ljl [l!| . Given 
the puncture motion frequency uj p (t) and the frequency 
of a comparable zero-eccentricity simulation u) c (t) (esti- 
mated by fitting a fourth-order polynomial in t through 
the numerical w p (t), as suggested in [3), the eccen- 
tricity is estimated by finding extrema in the function 
(u) p (t) — uj c (t))/(2uj c (t)). The uncertainty in the eccen- 
tricity estimate is about 2 x 10~ 4 27]. A simulation 
with initial D = 12M but using "quasi-circular orbit" 
parameters (as discussed further in Section TV C[) has an 
eccentricity of e w 0.008, i.e., five times larger than the 
eccentricity of the D12 simulation. 

For comparison between PN and numerical results, we 
must make clear what we mean by the individual black 
hole masses Mi and M2, and the total mass M. The 
mass of each black hole, Mi, is specified in terms of the 
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Arnowitt-Deser-Misner (ADM) mass at each puncture. 
This corresponds to the mass at the other asymptotically 
flat end, and has been found to equal numerically the 
apparent-horizon mass [37| , which for nonspinning black 
holes is related to the area of the horizon Ai by 

m > - M- m 

We assume that this mass is the same as the mass used 
in post-Newtonian formulas. Rather than try to quantify 
the accuracy of this assumption, we make the following 
argument, which we consider to be a more practical ap- 
proach. Our assumption is rigorously true only in the 
limit where the black holes are infinitely far apart and 
stationary. As such we consider any error in this assump- 
tion as part of the error due to starting the simulation at 
a finite separation. Since there is no invariant measure of 
quasi-local mass in general relativity, this error is present 
in some form in all numerical simulations. In practice one 
could rescale the total mass to optimize the phase and 
amplitude agreement with post-Newtonian calculations, 
but in the present work we retain the assumption that 
the horizon mass and PN masses can be equated. This 
provides an overall scale M = Mi + M2 for both numer- 
ical and post-Newtonian waveforms, and is crucial for 
comparison and matching. 

Let us discuss some other possible sources of error re- 
lated to the masses and spins of the black holes in our 
numerical simulations. 

The initial data contain "junk" radiation that quickly 
leaves the system. Some of this radiation may fall into 
the black holes and alter their masses. To estimate this 
effect, we refer to the initial-data studies of Cook and 
York [3a . |39| . who estimated the maximum radiation con- 
tent of single boosted Bowen-York black-hole initial-data 
sets (recall that a single boosted Schwarzschild black- 
hole spacetime will not contain any gravitational radia- 
tion). An estimate based on their data suggests that the 
spacetime of a Bowen-York black- hole with Pi /Mi « 0.17 
(which is the case for the D12 simulations) has a maxi- 
mum gravitational-wave energy content of 0.01% of the 
mass. In our simulations, the radiated energy from the 
junk radiation is at least 0.005% of the initial mass. 
Therefore we estimate that at most 0.005% of the mass 
fell back into the black hole. An error in our estimate 
of the total mass of 0.005% would lead to a phase er- 
ror in a 2000Af simulation of 0.1M. We calculate (See 
Section IIV|) a numerical uncertainty in the merger time 
of 0.4M, making any effect due to junk radiation falling 
into the black holes lower than our numerical uncertainty, 
and therefore not detectable at our level of accuracy. 

A further possible issue with the mass is that it may 
drift due to numerical error over the course of the sim- 
ulation. However, since we see clean sixth-order conver- 
gence in the time when the gravitational wave reaches 
a maximum, we expect that any mass drift either also 
converges away at sixth-order, or is well below the error 
due to other numerical effects. 



Finally, one may worry that the black holes pick up 
spin during their evolution. This effect has already been 
studied by Campanelli, et. al. [S(|. We do not at- 
tempt to measure this effect in our simulations, for the 
following reason: we are comparing numerical and PN 
waveforms of binaries that initially consist of nonspin- 
ning black holes. In the PN approach we use, the black 
holes remain nonspinning. Any spin that they acquire 
in full general relativity will therefore contribute to the 
disagreement between PN and NR waveforms. It is that 
difference in the waveforms that we are interested in mea- 
suring. More detailed investigation of the physical prop- 
erties of nonspinning binaries is beyond the scope of this 
study. 



III. POST-NEWTONIAN WAVEFORMS 

Binary inspiral waveforms can be constructed by a va- 
riety of means. We choose to compare our numerical 
waveforms with particular PN waveforms that are pro- 
posed for future searches for gravitational wave signals 
from black hole binary coalescence, namely the Taylor- 
expanded or EOB-resummed waveforms imp lemented in 
the LSC Applications Library (LAL) [1, 5J HI] . In par- 
ticular, we compare with the 3.5PN Taylor Tl waves, 
with a version of the code 1 that includes modifications 
to the flux coefficients given in the Erratum to [U . 
In the Taylor Tl approach ordinary differential equations 
are solved numerically to give the phase of the wave, and 
the amplitude is estimated by the quadrupolc contribu- 
tion, which is proportional to x = (Mw/2) 2 ' 3 , where ui 
is the gravitational- wave frequency, and lo/2 is assumed 
to be the orbital frequency of the binary. This treatment 
of the amplitude yields "restricted" PN waveforms. In 
Section [V] we also compare with a 2.5PN treatment of 
the amplitude [43[, which includes terms up to x 7 / 2 . 

To check the consistency of our comparison, we also 
compare with the "Taylor T3" PN approximant [H, 
[4f|, which consists of an analytic expression for the 
gravitational-wave phase as a function of the variable 
t = v(t — t c )/(5M), where t c is the "coalescence time" 
of the binary, M is the total mass, and 77 = M1M2/M 2 
is the symmetric mass ratio. The T3 approximant for 
the phase also contains a free phase constant, c/)q. The 
coalescence time t c and phase constant <po can be chosen 
to line up the phase and frequency of a T3 PN waveform 
with an NR waveform at a given time. We use the Tay- 
lorT3 approximant up to 3PN order, because the 3.5PN 
term contains an unphysical turning point long before the 
merger, which was already noted in the PN comparisons 
made in [HB1. 

The LAL code that we use, LALInspiralTest, pro- 



We used a version of LAL consistent with cvs version 1.25; earlier 
versions contain errors in the TaylorTl implementation. 
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duces h + and/or /i_ as a function of time. From this 
function we can compute the real part of ^ , 4 i (; = 2,m=2) 
by differentiating twice with respect to time. We choose 
^4,22 as the quantity to compare between NR and PN 
waveforms for two reasons: (1) we can compute the PN 
with arbitrarily small discretization error, and thus 
expect that its derivatives will be more accurate than 
computing by integration of the NR ^4,22! (2) in- 

tegration of ^4,22 requires estimating two constants of 
integration, which further complicates the procedure. In 
short, it should be equivalent to compare waveforms us- 
ing and ^4,22, and we choose ^4,22 because it is 
more straightforward. 

To generate a PN waveform we must choose the masses 
of the two bodies, and a range of frequencies that we 
want the waveform to cover. The masses are specified 
in units of solar mass. To produce the quantity 1^4,22 
that we wish to compare with numerical data, the time 
is rescaled to be in units of M by multiplying by the 
factor c/M, where we chose M = 2M Q = 2953.25 m 
(although the choice of masses is arbitrary) and the speed 
of light is c = 2.9979 x 10 8 m/s. The 3.5PN wave strain 
is then differentiated twice with respect to time, and the 
amplitude is scaled by the factor ^l6n/5 to give the 
coefficient of the I = 2, m = 2 mode. 

IV. NUMERICAL SIMULATIONS: ACCURACY 
AND CONSISTENCY 

In this section we describe our procedure for producing 
the most accurate waveform possible from our numerical 
data. This consists of taking waveforms calculated at five 
extraction radii from simulations performed at three res- 
olutions, and (1) Richardson extrapolating these wave- 
forms with respect to numerical resolution to produce 
accurate waveforms at each of the five extraction radii, 
and then (2) extrapolating with respect to extraction ra- 
dius to estimate the signal that would be calculated as 
R ex — > co. 

In [r| we described the use of sixth-order accurate spa- 
tial finite differencing in the bulk in BAM, introduced to 
increase the overall accuracy and in particular reduce the 
phase error in long evolutions. We found that Re(r^i)22, 
as directly computed by the code, was sixth-order con- 
vergent only up to about lOOAf before merger. However, 
if we separate the waveform into its amplitude and fre- 
quency as 

r* 4 = i4(#t))e*« (2) 

and examine separately <f>(t) and A(cj>), then the phase 
shows reasonably clean sixth-order convergence through- 
out the evolution (with a small "blip" around the merger 
time), and the amplitude computed as a function of the 
phase angle shows good convergence with far lower er- 
rors than when we consider simply A(t). This is because 
A(t) includes errors from the phase as well as the ampli- 
tude measurement; considering A(<j>) allows us to isolate 



the phase errors from the amplitude errors. With this 
phase/amplitude split we are able to perform Richardson 
extrapolation and reconstruct a more accurate waveform 
and calculate an error estimate. More details about the 
convergence properties of these simulations can be found 

in m. 

In order to be as clear as possible about this proce- 
dure, we will outline in detail the steps we followed to 
produce the D12 waveform that will form the basis of 
our comparison with PN waveforms. 

We perform three simulations with the grid configura- 
tion (following the notation in [l3|) x?j=2[5 x N : 5 x 2N : 
6], where N = 64, 72,80 for the D12 runs. The grid res- 
olutions on the finest inner box are Af/42.67, M/48 and 
M/53.33, and the resolutions on the coarsest outer levels 
are 12M, 10.67M and 9.6M, placing the outer bound- 
ary at about 775 M. The wave extraction is performed 
at resolutions 1.5Af, 1.33A/ and 1.2M. The grid setup 
is summarized in Table HI which also provides the grid 
details for D8, D9, D10 and Dll simulations. 

In each simulation, waves are extracted at radii R ex — 
40,50,60,80 and 90Af. Figure □ shows that the phase 
displays good sixth-order convergence over the course of 
the entire evolution. In order to disentangle the error 
in the phase from that in the amplitude, we now con- 
sider the amplitude as a function of phase, rather than 
time, A(4>), and show in Figure [5] that this function is 
also sixth-order convergent. For comparison Figure[3]also 
shows a convergence plot of the amplitude as a function 
of time, A(t), with no adjustments made for the phase. 
We see that A(t) is sixth-order convergent, but the errors 
are almost a factor of ten larger than they are for A((f>); 
this demonstrates the utility of considering A(<f>) instead 
ofA(t). 

The figures show the amplitude and phase from the 
waves extracted at R ex — 60Af, but similar properties 
are seen at all five extraction radii. Note that in these 
figures, and in all other relevant figures in this paper, 
the horizontal axis displays the time from the numerical 
code. For example, in Figure [T] the wave phase shown 
at t = 1000 M is the phase of the wave measured at the 
extraction sphere at R ex — 60M at code time t = 1000A/. 
In subsequent plots, when some time shifting has been 
applied, we indicate how this relates to the code time as 
displayed in any figures. 

Given the clean sixth-order convergence of A(jf) and 
4>(t), we apply Richardson extrapolation to A{<j>) and <fi(t) 
at each extraction radius. Since we have results at three 
resolutions, we are also able to compute an error estimate 
for the Richardson-extrapolated results. If a function in 
the continuum limit is /, and a numerical calculation of 
it, /, is sixth-order accurate, then we can write 

f = f + ai h 6 +a 2 h r + 0(h s ), (3) 

where h is the grid-spacing. With results at two resolu- 
tions, Richardson-extrapolation involves calculating the 
coefficient a\ and removing the sixth-order error to give 
a result that is seventh-order accurate. With results at 
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FIG. 1: Convergence of the phase <j>(t). Differences between 
simulations with N = 64, 72, 80 (see Tabled are scaled assum- 
ing sixth-order convergence. The convergence of the phase is 
shown as both a standard and a logarithmic plot, to demon- 
strate that good sixth-order convergence is seen throughout 
the simulation, except after merger, when there is a slight 
drop in convergence. In the logarithmic plot the solid and 
dashed lines are so close as to be almost indistinguishable. 



three resolutions, we may also calculate 02, and taking 
the difference between estimates of the true solution / 
using only a,\ or both a\ and 02, we can estimate the er- 
ror in the Richardson-extrapolated result. These errors 
are shown in Figure 2] for the portion of the simulation 
that will be compared with PN waveforms. We see that 
for t > 500M the uncertainty in <j>(i) is less than 0.01 
radians, and the uncertainty in A((f>) is less than 0.5 per- 
cent. At earlier times the uncertainties grow by up to a 
factor of ten, due to noise in the data. 

We now have amplitude and phase functions A((f>) and 
(j)(t) for each of the five extraction radii, and wish to 
extrapolate to R ex — > 00. 

We first deal with A(4>). Since we are looking at the 
amplitude as a function of phase, rather than time, the 
amplitudes measured at each extraction radius are al- 
ready in phase; there is no need to "line them up" , as 
would be necessary if we looked at A(t). We find that 
the value of the five amplitude functions is approximated 
well by a quadratic function in extraction radius, i.e., 



A(4>,R ax )=A co (4>) + 



R 2 



o 



1 



(4) 



In other words, the wave amplitude falls off as the square 
of the extraction radius. A simple curve fit (performed at 
each phase <f) allows us to construct A^ty). Including 
the next fall-off term, 1/R^ X , allows us to also estimate 
the uncertainty in the extrapolation, analogous to the 



FIG. 2: Convergence of the amplitude A(<f>). Differences be- 
tween simulations with TV = 64, 72, 80 (see TableHJ) are scaled 
assuming sixth-order convergence. The a;-axis shows <^/(47r), 
which gives a rough estimate of the number of orbits the sys- 
tem has completed (at least before merger). The phase <f> is 
chosen to be zero at t — 0. The convergence of the amplitude 
is shown in terms of relative (percentage) errors, to allow eas- 
ier comparison with later results. A vertical line indicates 
the point at which we end our PN comparison in Section IVl 
The lower plot zooms into the region that will be used for PN 
comparison. 
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FIG. 3: Same as Figure [21 but using A(t) instead of A(<f)). 
We see that the errors are far larger than for A(cj>); the maxi- 
mum error is now around 60%, while it was only 8% when we 
considered A(cj>). 



method of error estimate in the Richardson extrapolation 
of the discretization error. Note that although one would 
expect the error to fall off as 1/ R ex , our results suggest 
that the quadratic fall-off dominates; this has also been 
observed in simulations of a particle orbiting a Kerr black 
hole [46[ . The quadratic fall-off in the amplitude error is 
demonstrated in Figure [5] The resulting relative error 
estimate as a function of </> is shown in Figure [51 and as 
a result of this plot we estimate the uncertainty in A((j>) 
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FIG. 4: Error in the Richardson-extrapolated functions <f>(t) 
and A{4>). For the range of the simulations that will be com- 
pared with PN waveforms, the uncertainty in <f>(t) is below 
0.01 radians at most times, and the uncertainty in the ampli- 
tude is less than 0.5%. At earlier times (t < 500M, which are 
also nominally included in the PN comparison), these plots 
are dominated by noise and the uncertainty grows by a factor 
of ten. 
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FIG. 5: The wave amplitude A as a function of extraction 
radius R ex , at 4> = 87r, which corresponds to t w 715M for the 
wave extracted at Rex = 90M. The solid line shows a curve 
fit of the form ([4]). The dashed line shows a curve fit with an 
extra 1/R^x term. The horizontal solid and dashed lines show 
the corresponding R ex — > oo limits of the two curve fits; our 
uncertainty estimate in the extrapolation of the amplitude 
comprises the difference of these two values. 



due to extrapolation to R ex — > oo as about 2%. This 
dominates the uncertainty from Richardson extrapola- 
tion (< 0.5%), so we also estimate the overall uncertainty 
in A(4>) as about 2%. 

We now turn to the phase, 4>{t). To a first approxima- 
tion we expect that the difference in the phase measured 
at two extraction spheres will be a constant. However, 
the proper distance between each extraction sphere may 
drift due to gauge effects. We have already seen in evolu- 
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FIG. 6: Error in the R e x — > oo extrapolated function A(<f>). 
For the range of the simulations that will be compared with 
PN waveforms, the uncertainty in the amplitude is less than 
2%. 



tions of the Schwarzschild spacetime that the coordinate 
location of the horizon drifts depending on the value of 
the r\ parameter in the T-driver shift evolution equation 
(see Figure 4 in (HI), and effects related to r\ have also 
been observed and studied in [I?], |H| ; and it is quite pos- 
sible that there are other gauge effects that we are not 
aware of. 

We have attempted to extrapolate the phase to R ex — > 
oo by lining up the phase at a given time, and then ob- 
serving, at other times, the deviations in the phase at dif- 
ferent extraction radii. These deviations decrease as R ex 
increases, and the fall-off can be reasonably well modeled 
by a polynomial in 1/ R ex , and far better by a polynomial 
in I /Rex and \/R\ x . However, we do not find the limit as 
R ex — > oo to be very robust — the results vary depending 
on the choice of the time when the phases are lined up. 
(Obvious choices for this time are when the gravitational 
wave amplitude reaches a maximum, near merger, or the 
time at which the GW frequency Moj equals one of the 
matching values that will be used in our PN comparison 
below.) As such, we do not extrapolate the phase. We 
instead use the phase at the largest extraction radius, 
Rex = 90M (which we expect to be the most accurate) 
and use the phase extrapolation procedure to estimate 
the uncertainty in the phase, which we give as 0.25 radi- 
ans. 

An alternative indication of the accumulated phase er- 
ror of the numerical simulations is given by the time when 
the amplitude of the gravitational-wave signal reaches a 
maximum. This time is also seen to be sixth-order con- 
vergent, and a similar Richardson-extrapolation error es- 
timate as performed above gives an uncertainty of 0.4M 
in the "length" of the simulation. 

We are now able to construct a final waveform, 



i?e(r^ 4 ,2 2 (<)) = A c 
Jm(r* 4 ,22(*)) = A c 



,(0 9O (<))cos(0 9O (<)-5) (5) 
,(09o(t))si]i(09o(t) - 5), (6) 



where 6 is an arbitrary phase shift, which we will ap- 
ply later when comparing with PN waveforms. The un- 
certainty in the wave amplitude is about 2%, and the 
accumulated phase error over the time range we will 
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consider is about 0.25 radians. The time-shifting pro- 
cess described earlier means that the extrapolated wave- 
form is measured at an effective extraction sphere with 
Rex — 90M, i.e., our extrapolated waveform gives the 
wave amplitude that would be measured at infinity, but 
at a time roughly 90Af after the wave was emitted from 
the binary system. Since we do not make direct com- 
parisons between quantities calculated from the gravita- 
tional waves and quantities calculated from the puncture 
motion, it is not necessary to know this "wave travel- 
time" precisely. Although not used here, one could esti- 
mate this time using the method suggested in [4!| . 

The same procedure is applied to the simulations D10 
and Dll. With a suitable time-shift applied so that the 
wave amplitude maximum occurs at the same time, the 
extrapolated waveforms from the three simulations are 
shown in Figure[71 The waveforms lie almost perfectly on 
top of each other, except in the last cycle before merger. 
It it at this time that we see a "glitch" in the clean con- 
vergence of the phase <fr(t). However, for comparison with 
3.5PN waveforms we will only be interested in the wave- 
form before t w 1770 M . In order to quantify the level of 
agreement between the D10, Dll and D12 waveforms, we 
also show in the lower panel of Figure [7] the accumulated 
phase error between t = 1200M and t = 1800Af (where 
t is the code time of the D12 simulation). We see that 
the phase errors average to below 0.03 radians. 



V. COMPARISON 

Given the post-Newtonian (PN) and numerical- 
relativity (NR) waveforms discussed in Sections IIIII and 
IIV1 we are now in a position to compare them. We com- 
pare NR waveforms with a 3.5PN TaylorTl waveform 
that was terminated at a gravitational-wave frequency 
Moj = 0.120, but we will only use it up to a cutoff fre- 
quency of Mu>o =0.1; since the growth in phase error in 
the 3.5PN waveform becomes dramatic at late times (see 
for example Figure 17 in [2lj]), the smaller the choice of 
cutoff frequency the better. Figure [8] shows the numeri- 
cal D12 r^ , 4 ! 22 overlaid with the 3.5PN TaylorTl version 
computed from output from the LAL code. The figure 
starts at t = 340Af , after the binary has completed one 
orbit; this allows time for noise due to the junk radiation 
in the initial data to leave the signal. The agreement 
between the PN and numerical waveforms appears to be 
excellent. A similar plot (for h+) is shown in Figure 1 of 

m. 

The NR and PN waveforms shown in Figure [5] were 
"lined up" by first identifying the time at which both 
waveforms had a given frequency Mwq. An appropriate 
phase shift S was then applied to the numerical waveform 
to line up the PN and NR phases. The choice of Mujq 
can have a dramatic effect on the quality of the phase 
agreement between the PN and NR waveforms. Figure [5] 
was produced by matching at the beginning of the com- 
parison region, at Mcu = 0.0455, which gives a far better 
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FIG. 7: Final waveforms from the D10, Dll and D12 runs, 
produced by the method discussed in the text and shifted in 
time so that their amplitude maxima occur at the same time. 
The three lines are not individually labeled; the main point is 
that their differences are almost indistinguishable, except in 
the last cycle before merger; see text. The phase disagreement 
with the D12 simulation is shown in the lower panel. 

phase match, as we will discuss below. 

We will now discuss this subtle feature of the matching 
process in more detail, before we make any conclusions 
about the agreement between NR and 3.5PN TaylorTl 
waveforms. 



A. Phase and frequency 

The wave frequency Mui calculated from the NR waves 
is typically very noisy at early times, but becomes much 
smoother near merger, when the value is higher. To al- 
low a matching at any time in the window of compari- 
son, we fit a polynomial in time through the numerical 
frequency to produce a smoother function. The curve fit 
is based on the form of the frequency evolution in the 
TaylorT3 approach, i.e., a polynomial in (t — t c ), where 
t c is a crude estimate of the merger time (its specific 
value does not strongly affect the accuracy of the fit; we 
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FIG. 8: Numerical (solid line) and TaylorTl 3.5PN (dashed 
line) waveforms for equal-mass inspiral. 



used t c — 1927M), and the powers of [t — t c ) that are 
included are {-3/8,-5/8,-3/4,-7/8,-1,-9/8}. The 
use of a curve fit introduces yet another source of error 
in our numerical phase, particularly at early times, which 
is difficult to assess. However, the analyses below were 
repeated with different fitting functions (by keeping or 
removing the last term in the fit, or varying t c ), and all 
changes in the phase results were below the stated nu- 
merical phase uncertainty of 0.25 radians. Nonetheless, 
we tend to consider any matching done at late times to 
be more reliable than that done at early times. 

On the other hand, we expect the PN phase to be 
most accurate at early times — in principle, we should be 
able to obtain arbitrary accuracy in the post-Newtonian 
expressions by going to sufficiently early times. For that 
reason we first choose to line up the frequencies at t = 
347AM in code time (recalling that this is the time when 
the wave reaches the extraction sphere at R ex = 90M), 
when Mu) = 0.0455. We are then free to make a constant 
phase shift S to align the phase of the waves; again aligned 
at t = 347. AM with 6 = 1.3677T. The agreement between 
the NR and 3.5PN wave frequencies as a function of time 
is shown in Figure [9l As can be seen in the lower panel 
of Figure [51 the PN and NR frequencies remain close up 
to around t = 1000M, and then drift apart and finally 
diverge. Also shown (with a dashed curve) is the result 
of matching at the end of the comparison region (at t — 
1772M, with Mu = 0.1, S = 1.067tt). 

The corresponding results for the phase disagreement 
are shown in Figure [TOJ Also shown is the phase dis- 
agreement between the NR waveform and a waveform 
produced using the TaylorT3 approximant. In order to 
line up the phase and frequency of the T3 waveform, we 
choose an appropriate coalescence time t c and phase con- 
stant <j)Q. 

Figure [10] demonstrates that the different choices of 
matching frequency can give entirely different impres- 
sions of the relative merits of the Tl and T3 approxi- 
mants: when the waves are matched at t = 1772M, the 
accumulated phase disagreement between the T3 approx- 
imant and numerical results is about 0.1 radians. When 
the matching is done at t = 347AM, the accumulated 
T3/NR phase disagreement is almost 1 radian. In both 
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FIG. 9: The wave frequency u as a function of time for 
the D12 numerical-relativity and TaylorTl 3.5PN waveforms. 
The frequencies agree at t = 347AM, when Mlo = 0.0455 
The percentage disagreement between the two is shown in 
the lower panel. Also shown is the frequency disagreement 
when matching is done at Mw = 0.1, t = 1772M. 

cases the Tl/NR disagreement is comparable, although 
this is purely an accident of the matching frequencies that 
were chosen. It should be clear from the lower panel of 
Figure [TU] that if we cut off the comparison at t = 1000A/, 
the Tl/NR accumulated phase error will be very small. 
Similarly, for matching purposes, one could optimise the 
matching time to give the smallest phase error — for the 
Tl waveforms, we can for example match at Mlu s=s 0.075 
and achieve a phase agreement within numerical uncer- 
tainty. 

We repeat that for the purposes of comparing PN and 
NR phases, the match at Mcu =0.1, when the numerical 
data is relatively free of noise, is the most trustworthy. 
The matches at earlier times are less accurate and mainly 
serve to illustrate the general trend in the disagreement 
between PN and NR phases: the frequency disagreement 
changes sign (as shown in Figure 0, and, depending on 
the approximant used and the chosen matching time and 
frequency, the phase disagreement may behave as in the 
T3/NR curve in the top panel of Figure [TO] or the Tl/NR 
curve in in the lower panel, and exhibit a local maximum, 
which allows us to optimize the phase disagreement. 

We may produce yet another picture of how Tl and 
T3 behave by plotting the phase disagreement versus the 
wave frequency Mu, as done in [15j. This is shown in 
Figure QTJ which now suggests that T3 behaves far better 
than Tl. 

What are we to conclude, then, about the phase agree- 
ment between NR and Tl or T3 PN waveforms? Due to 
a turning point in the evolution of the frequency dia- 
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FIG. 10: The disagreement in the phase between NR wave- 
forms and PN waveforms constructed with the TaylorTl and 
TaylorT3 approximants. In the upper plot the phase and fre- 
quency are matched at t = 1772M, Muo = 0.1. In the lower 
plot they are matched at t — 347AM, Mui = 0.0455. We see 
that the relative merits of the two approximants can appear 
quite different depending on the matching time. 
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FIG. 11: The disagreement in the phase between NR wave- 
forms and PN waveforms constructed with the TaylorTl and 
Taylor T3 approximants, but now shown as a function of GW 
frequency Mui. 



greement, we are left with a great deal of freedom about 
how to match the frequency and phase. We find, in the 
frequency range that we consider, that the minimum ac- 
cumulatccd phase disagreement that we can achieve is 
about 0.2 (or 0.15) radians using either the Tl (or T3) 
approximants (see Figure fTO]) . By contrast, the maxi- 
mum disagreement between the NR and PN phases over 
the comparison region is about 1 radian, although since 
this results from a matching at early times, and the phase 
disagreement is diverging at the end of the comparison 
region, this value has a large uncertainty. 

In a matching between NR and PN waveforms (as per- 
formed in, for example, @, we naturally choose to 
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FIG. 12: NR and restricted PN amplitudes of r* 4 ,: 



match in such a way that the phase disagreement is min- 
imized. We could easily have found that the PN phase 
evolution disagreed so badly with the NR phase evolu- 
tion that it was not possible to achieve an accumulated 
phase disagreement of less than, for example, 1 radian. 
However, the minimum accumulated phase disagreement 
that we can achieve is about 0.2 radians, which is also 
within the phase uncertainty of the numerical waveforms. 

We therefore conclude that we can match the phase 
within the numerical uncertainty over the frequency 
range we have considered (Mui — 0.0455 up to Mui = 
0.1), and that the accumulated PN and NR phase dis- 
agreement has an upper bound of roughly 1 radian. We 
expect that matching at even earlier times (using longer 
simulations) would make the matching clearer, although 
this will also require more accurate simulations and larger 
radiation extraction radii to resolve the lower-frequency, 
lower amplitude waves. 



B. Amplitude 

We now turn to the amplitude. 

Figure [T2l shows the amplitude of r\&4 i 22 from NR and 
restricted PN waves, plotted as a function of GW fre- 
quency Mui, so that the choice of PN approximant does 
not affect the result. The amplitude of the restricted 
3.5PN wave is larger than that for the NR wave. Fig- 
ure [T2] shows the percentage disagreement between the 
restricted PN and NR wave amplitudes over the same 
frequency range. The disagreement is of the order of 
6%. Since the uncertainty in the NR wave amplitude is 
below 2%, at least for Mui > 0.05, we cannot ascribe 
this disagreement entirely to numerical error. If we as- 
sume that the NR wave more closely models the cor- 
rect physics of the binary system, then the restricted PN 
(quadrupole) amplitude over-estimates the amplitude by 
between 4 and 8% in this frequency range. 

So far we have compared our NR waveforms with 
restricted 3.5PN waveforms, meaning that the ampli- 
tude in the gravitational-wave strain is proportional to 
x = (Mw/2) 2//3 . (The factor of two signifies that x deals 
with the frequency of the black holes' motion, not the 
frequency of the waves; the two frequencies are assumed 
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FIG. 13: Percentage disagreement between the NR wave 
amplitude and that of the PN waves with the restricted 
(quadrupole) and 2.5PN amplitude treatments. At low fre- 
quencies the 2.5PN amplitude agrees with the NR amplitude 
within numerical uncertainty. 
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FIG. 14: The percentage difference in amplitude between the 
D12 simulation and the D9, D10 and Dll simulations. The 
D9, D10 and Dll amplitudes show small oscillations around 
the D12 value, but recall that the disagreement between the 
D12 and restricted 3.5PN amplitudes was 6%. 



to be related by a factor of two.) If we move beyond 
restricted waveforms, and model the amplitude up to 
2.5PN order (i.e., with terms up to x 7 ^ 2 ) [l5|,[43|, we find 
greater disagreement at higher frequencies, but at low 
frequencies the 2.5PN amplitude shows better agreement 
with the NR amplitude. The 2.5PN amplitude disagree- 
ment at Mui = 0.0455 is between 1% and 5%; the PN and 
NR amplitudes now agree within numerical uncertainty. 
This is also shown in Figure [T3J 

As we have said, the amplitude disagreement between 
the NR and restricted PN amplitudes is roughly constant 
over the frequency range Mu> — 0.05 to Mu> — 0.1. This 
suggests that if we are content with these levels of error 
when matching numerical and PN waveforms, the large 
number of cycles in the D12 simulation is not necessary. 
A combined PN-NR waveform could be produced by ap- 
plying a scale factor, as is done using different approaches 
in [2J, [2f| and [2l[ , and clearly only a few cycles shared 
by the NR and 3.5PN waveforms are needed to deter- 
mine the scale factor. We may now ask: can we get away 
with a numerical simulation that starts at, for example, 
D = 9M, and yields a waveform that (neglecting the first 
orbit) shares four cycles with the 3.5PN wave? 

Figure [T4l shows the relative disagreement in amplitude 
between the D12 simulation and the D9, D10 and Dll 
simulations. There are small oscillations around the D12 
values, but these are smaller than the average amplitude 
disagreement between the NR and restricted 3.5PN wave 
amplitudes, and we expect that it will be possible to cal- 
culate a suitable scale factor for matching the NR and 
3.5PN waves. We conclude then that simulations start- 
ing as close as D — 9M and simulating about 4.5 orbits 
should be enough to match to restricted 3.5PN wave- 
forms for many applications. To make this clearer: any 
GW data analysis application that requires an amplitude 
accuracy of at most 5% up to the last four cycles, and 
an amplitude accuracy of better than 2% from that point 
through merger and ringdown, will require only short (4.5 
orbit) numerical simulations to match to PN waveforms. 

This result is attractive from a computational point of 
view. The D9 simulation ran in 750 CPU hours (two and 



a half days of wall clock time on 12 processors), while the 
highest resolution D12 simulation required 10,500 CPU 
hours (18 days on 24 processors). When producing many 
waveforms for use in gravitational-wave data analysis, we 
would much rather only have to perform the two-and-a- 
half-day simulations. 

Of course, in the case of equal-mass binary inspiral, 
we have already presented waveforms that cover far more 
than four cycles before merger. The important question 
is whether similarly short simulations will be adequate 
beyond the equal-mass nonspinning case, and that will 
be the subject of future work. 

C. Comparison with eccentric waveforms 

The numerical simulations discussed in the previous 
sections modeled equal-mass inspiral with negligible ec- 
centricity, starting from the initial parameters introduced 
in [27[. The eccentricity of the D12 simulation is esti- 
mated as e < 0.0016. In contrast, one could use standard 
"quasi-circular orbit" parameters (i.e., parameters calcu- 
lated with the assumption that f = 0), which lead to 
inspiral with a small but noticeable eccentricity. We now 
consider a set of simulations with the same parameters as 
the D12 runs, but using initial parameters calculated us- 
ing the 2PN-accurate expression used in [5(| (and based 
on the results in [5l|); we denote this simulation "QC12". 
We apply the same extrapolation procedure as described 
in Section IIVI to produce the final waveform that we an- 
alyze. 

Figure [T5] shows the same comparison with the Tay- 
lorTl 3.5PN wave phase as in the upper panel of Fig- 
ure[T0l but now displaying results from both the D12 and 
QC12 simulations. The accumulated phase disagreement 
for the QC12 simulation is larger. The disagreement with 
the 3.5PN phase also shows oscillations that are presum- 
ably due to eccentricity. A similar effect can be seen 
in Figure 1161 which shows the percentage disagreement 
in wave amplitude. The amplitudes are now shown as 
functions of time; if we use Muj as before, the eccen- 
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FIG. 15: The same quantities as in Figure [9] this time com- 
paring both the D12 and QC12 wave phase with that of the 
TaylorTl 3.5PN waves. The disagreement between the 3.5PN 
and QC phases displays clear oscillations, presumably due to 
the eccentricity in the QC12 simulation. 



12 
10 

£ 8 
I 6 
XJ 4 



Time (M) 

FIG. 16: Percentage disagreement between restricted 3.5PN 
and NR wave amplitudes, for both D12 and QC12 simula- 
tions. The disagreement between the QC12 and 3.5PN wave 
amplitudes is clearly dominated by eccentricity. The low- 
eccentricity D12 simulation is necessary to identify the error 
in the restricted 3.5PN wave amplitude. 
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D. Comparison of the black hole coordinate motion 

To initialize our numerical simulations, we have set the 
initial momenta of the black holes to values we have ob- 
tained from a post- Newtonian inspiral calculation as de- 
scribed in 27[ ■ The inspiral calculation starts at an initial 
separation of D — 40M with momenta given by the 3PN- 
accurate quasicircular-orbit formula given in [131 ] . When 
the inspiral reaches the separation D = 12M, the mo- 
menta are read off from the solution and given the values 
shown in Table [II] 

In this section we compare a full GR simulation that 
uses those parameters with continuing the PN inspiral 
from D = 12M. 

The coordinate separation of the black-hole punc- 
tures was chosen as the coordinate separation of the 
post-Newtonian inspiral, which we have computed in 
ADMTT-coordinates. This is motivated by the fact that 
the PN solution in the ADMTT gauge for a two-body sys- 
tem agrees with our Bowen-York puncture initial data up 
to 2PN order (see, for example, the explicit solutions in 
Appendix A of [52|). It is therefore interesting to know 
when the use of the ADMTT gauge breaks down in our 
evolutions. An indirect check is straightforward: we com- 
pare the PN and full NR puncture separation, as seen in 
Figures [IT] and [18] Using the D12 simulation, we find 
that both the separation and orbital phase agree very 
well from D = 11 M up to D = 8M, or from t = 300M 
(the time to complete the first orbit) to t — 1500A7. Put 
another way, the PN and full NR coordinate separation 
agrees until about 3 orbits before merger. 



VI. DISCUSSION 



tricity effects are not visible. The low-eccentricity D12 
waveform has been matched with the PN waveform at 
Mu> = 0.0455, and the QC12 waveform is matched with 
the D12 waveform so that their amplitude maxima occur 
at the same time. The amplitude disagreement between 
the D12 simulation and the 2.5PN amplitude is slightly 
different to that shown in Figure[l4] this is due to param- 
eterizing the amplitude with time instead of frequency — 
the PN/NR frequency disagreement means that there is 
not a 1-1 relationship between the two plots. However, 
the results are consistent within the 2% uncertainty in 
the numerical waveform amplitude. 

The disagreement in amplitude between the 3.5PN and 
QC12 results oscillates between 2% and 10% at early 
times. From the QC12 simulation alone, we may guess 
that the error in the 3.5PN wave amplitude is the av- 
erage of this curve, i.e., around 6%, but may also guess 
that the disagreement might go away if the eccentricity 
were removed. The D12 simulation, which displays far 
less eccentricity, confirms the first guess: there is strong 
numerical evidence that the restricted 3.5PN wave am- 
plitude really does disagree with fully general-relativistic 
results by about 6%. 



We have simulated nine orbits, merger and ringdown 
of an equal-mass binary, and extracted waveforms of 
sufficient accuracy to make a detailed comparison with 
post-Newtonian (PN) waveforms. The uncertainties in 
the numerical waveforms are dominated by the close ex- 
traction radii, and not finite-difference errors. The PN 
waveforms that we focused on were those generated by 
the TaylorTl 3.5PN procedure implemented in the LSC 
Applications Library (LAL), which is a candidate for 
use in gravitational-wave searches in detector data; we 
also compared with the TaylorT3 approximant. We find 
that the phase of the TaylorTl 3.5PN waveform can be 
matched to agree with the numerical phase to within nu- 
merical uncertainties, and the upper bound of the accu- 
mulated phase disagreement is on the order of 1 radian. 
The restricted PN amplitude overestimates the numer- 
ical value (6 ± 2)%. We have found that the ratio of 
the restricted PN and NR wave amplitudes is roughly 
constant over the course of the evolution, and therefore 
an equally good matching between PN and NR waves 
should be possible with far less numerical cycles. In par- 
ticular, we performed a simulation that completes only 
4.5 orbits before merger, and expect that this could be 
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FIG. 17: Orbital coordinate motion of the D12 numerical 
relativity evolution compared with a PN evolution with the 
same initial parameters. In both panels the PN evolution is 
drawn as a dashed line. Top panel: the separation of the black 
holes (the puncture position in the full NR case). Bottom 
panel: the coordinate angular velocity. 
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FIG. 18: Orbital coordinate motion of the D12 numerical 
relativity evolution compared with a PN evolution with the 
same initial parameters. Dashed line: (u>nr — lupn) /upn , full 
line: (D NR - D PN )/D PN 



matched to PN waveforms by a procedure like that dis- 
cussed in [23, [HI or [2l[ just as well as a simulation that 
models many more cycles. We therefore conclude that, 
with the level of numerical accuracy that we can achieve, 
only about 4.5 orbits need be simulated for a PN/NR 
matching of the same accuracy. Whether these relatively 
modest requirements for numerical waveforms carry over 
to the cases of unequal-mass and nonspinning binaries 
will be the subject of future work. 

For gravitational-wave detection we expect that such 
hybrid waveforms will be acceptable. However, for pa- 
rameter estimation the issue of the discrepancy between 
the amplitude of PN and NR waveforms may have to be 
addressed. Modeling the amplitude at 2.5PN order gives 



agreement within numerical error between PN and NR 
waves up to about 11 cycles before merger; at present we 
suggest that the best matching can be done with > 11 
cycles (5.5 orbits) of numerical simulation. The cases 
where the current level of phase and amplitude accuracy 
are expected to be adequate for various data-analysis ap- 
plications will also be explored in future work. 

Comparing with evolutions of the PN equations of mo- 
tion in the ADMTT gauge, we find that the orbital mo- 
tion seen in the numerical evolutions agrees extremely 
well up to a coordinate separation of about D = 8M. 
This surprising agreement not only suggests that the PN 
dynamics accurately models the full physical dynamics 
up to about three orbits before merger, but that the nu- 
merical gauge remains close to the ADMTT gauge up to 
that time. In addition, the gauge dynamics and emission 
of junk radiation at the beginning of the simulation do 
not noticeably change either the dynamics or the gauge; 
after about one orbit the NR dynamics matches up again 
with the ADMTT PN dynamics. 

Note: While this article was undergoing peer review, 
the Caltech/Cornell group completed a detailed PN/NR 
comparison that covers 30 gravitational-wave cycles (15 
orbits) before merger with high numerical accuracy in 
their numerical waveforms [53l |. Where comparable, their 
results confirm those in this paper; a comparison between 
our results and theirs' is provided in their paper. 
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